Avalanche dynamics in fluid imbibition near the depinning transition 
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We study avalanche dynamics and local activity of forced-flow imbibition fronts in disordered media. We 
focus on the front dynamics as the mean velocity v of the interface is decreased and the pinning state is ap- 
proached. Scaling arguments allow us to obtain the statistics of avalanche sizes and durations, which become 
power-law distributed due to the existence of a critical point at v = 0. Results are compared with phase-field 
numerical simulations. 
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In many physical systems, the response to a slow external 
driving usually involves avalanches or bursts. Different ex- 
amples are found in fracture cracks [1], granular material 10], 
earthquakes ft^, or during imbibition of fluids in porous me- 
dia Jj], among others. A particularly interesting problem in 
this context is the dynamics of fronts during imbibition of flu- 
ids in porous media [4], with its many engineering applica- 
tions in fluidics and oil-recovery technology. 

Imbibition in disordered media occurs when a viscous fluid, 
which wets the medium preferentially, displaces a less viscous 
fluid (typically air) and therefore, at relatively low injection 
rates, stable fronts separating the two phases are formed (for 
a recent review on imbibition see Ref. [51]). In the case of 
forced-flow imbibition, the spatially averaged velocity of the 
liquid-air interface v is kept constant by means of a constant 
injection rate. Then, for relatively low velocities, the invad- 
ing fluid advances in the form of spatially localized events or 
avalanches, as occurs in other disordered systems. 

Avalanche dynamics in imbibition is expected to be respon- 
sible for the front velocity fluctuations. Rost et al. |4J] have 
recently shown that in the case of imbibition, which is a lo- 
cally conserved process, velocity fluctuations are controlled 
by a length scale £ x arising from fluid conservation. This 
characteristic length scale introduces a natural cutoff in the 
distribution of avalanche sizes and durations, which leads to 
non-critical avalanche distributions and is ultimately respon- 
sible for the lack of correlated fluctuations at large distances. 
The scaling behavior of the velocity fluctuations can then be 
derived making use of the central limit theorem J4J] . Forced- 
flow imbibition described by avalanches with a fixed cutoff 
size has been experimentally observed in the recent work by 
Planet et al. [6], by means of analyzing the global velocity 
time series v(t). 

In this paper we study the statistics of local avalanches of 
activity in forced-flow imbibition in disordered media. We 
analyze the mesoscopic behavior of the interface by monitor- 
ing locally active sites, i.e. those sites that are moving at a 
given time, which define the actual avalanche taking place in 
the system (see Fig. Q]). We show that avalanche sizes and 
durations become power-law distributed for low enough in- 
jection rates due to the existence of a critical point at v = 0. 
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FIG. 1: (Color online) (a) Typical avalanches of the front h(x,t) 
from our phase-field simulations for v — vo/40 and clip c t h = 3. (b) 
Spatio-temporal activity is characterized by the size s, lateral extent 
I. and duration T of avalanches. 



The singularity appearing as v — > affects the value of the 
critical exponents that characterize the front dynamics and 
morphology- namely, the avalanche exponents, as well as the 
roughness exponents. This leads to effective exponents even 
for finite velocities. We obtain scaling relations connecting 
the roughness exponents and the avalanche exponents. Our 
scaling theory is compared with numerical results of a phase- 
field model for imbibition. 

Only in the last few years it has been possible to achieve a 
satisfactory theoretical understanding of imbibition, based on 
a detailed description of the physical forces that play a rele- 
vant role at different spatial scales |5J]. Surface tension a tends 
to flatten the front at short length scales, while quenched dis- 
order in both, capillary p c (r) and permeability K(r), makes 
the front to roughen and fluctuate around its average position. 
Both disorders operate at very different length scales, sepa- 
rated by a crossover length that depends on the velocity as 
£,k ~ 1/u Q5l LZD - We are concerned here with the interesting 
case of a slowly advancing front, i.e. the capillary dominated 
regime, where the permeability disorder is irrelevant and can 
be considered to be a constant K (r) ~ Kq. Different theo- 
retical approaches have arrived at the conclusion that, 
for small deviations around the mean interface position, the 
dynamical evolution for the liquid-air interface is described in 
Fourier space as 

d t h k = -<rK \k\k 2 h k - v\k\h k + K \k\i) k , (1) 

where fjk(h) are the Fourier components of the capillary dis- 
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order at some coarse-grained scale. From Eq. ([T]) one can eas- 
ily see that there exists a crossover length 
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such that interface fluctuations are uncorrelated above this 
typical scale. Indeed, several numerical studies Ji, lol 12] 
have shown that the interface is asymptotically flat on length 
scales larger than £ x , introducing then a natural cutoff in the 
system. For capillary-induced fluctuations we have £ x <C 
so that the permeability disorder can be ignored. 

Avalanche statistics.- In order to monitor local avalanches 
of forward movements we proceed as follows. First, we define 
the active sites on the interface as those where the local veloc- 
ity v(x, t) = dth(x, t) takes values above some fixed thresh- 
old, v(x, t) > cthv, where cth is some arbitrary constant and 
V is the spatially averaged global velocity. An avalanche is 
defined as a connected cluster of active [v(x, t) > c^v] sites 
surrounded by non-active [v(x, t) < Cthv] sites (see Fig. |TJ. 

Avalanches exhibit a typical size (volume) (s(£)) for an 
event of lateral spatial extent t. For a given front velocity 
v we expect the average avalanche size to scale with the lat- 
eral extent up to the cutoff length scale, (s(£)) ~ £ D for 
I <C £ x ~ v~ x l 2 , where D is the avalanche dimension ex- 
ponent that can be easily related with the local roughness ex- 
ponent a\ oc via the local width of the interface fluctuations 
w{£). One has (s(£)) ~ l d w{£) ~ £ d £ ai °% and D = d + ai oc 
in d + 1 dimensions. In particular, for d = 1, one observes 
that ai oc = 1 ll8l- [l0i[l3ll and then D = 2. We also expect to 
observe a scaling relation (s(T)) ~ T 1 between an avalanche 
of duration T and its size below a certain time cutoff T x . 

In order to study the statistics of the avalanche dynamics, 
we first calculate the probability densities V(s) and P(T) for 
having avalanches of size s and duration T, respectively. Due 
to the existence of the intrinsic crossover length in the imbi- 
bition problem the several avalanche probability distributions 
are not generically expected to be critical, but exponentially 
decaying functions: 



P(q) ~ q Te exp-{g/g x ), 



(3) 



where P(g) is the marginal probability density function 
(PDF); the index g denotes the size s, lateral extent I, or dura- 
tion T of avalanches, and r e is an exponent. The distribution 
cutoff depends explicitly on v and, in particular, the maximum 
avalanche size is given by 
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(4) 



This cutoff diverges (£ x — > oo) as the control parameter TJ — » 
0, which renders critical avalanches expanding over the whole 
system. This divergence is very strong, ,s x ~ (v)^ 1 , already 
in d = 1, which -in turn- is expected to be reflected in long- 
tailed avalanche distributions even for finite values of v (see 
Fig.0. 

Let us now consider the joint probability V(s,£,T) for 
having an avalanche of size s, extent £, and duration T. 
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FIG. 2: (Color online) (a) Distribution of avalanche size for dif- 
ferent velocities. The dashed curves correspond to a data fit to 
V(s) = as~ Ta exp [— s/s x («)]• The solid line fits the power-law 
regime for the smallest velocity. The inset shows the cutoff s x of 
the avalanche size distribution for each velocity with a guide-to-eye 
line of slope — 1. (b) Distribution of avalanche durations. The inset 
shows a comparison for different choices of the arbitrary threshold 
c t h from 3 to 20 in the case of v = vq /40 = 5x 10~ 5 a. u. 



At the critical point v = scale-invariant behavior implies 
V(s, I, T) = b a V{b D s, bl, b z T) for any scaling factor b > 0. 
Here, z corresponds to the interface dynamic exponent. In- 
tegrating over two of the arguments one obtains the marginal 
PDFs and the scaling relations 

t t = H (t s - 1), and t £ = 1 + z(t t - 1), (5) 

z 

that must be satisfied in the case of scale-invariant avalanche 
dynamics. They connect the avalanche activity exponents 
(tt, t s , T() with the dynamics of the front (z). Note that in 
this limit long-range interface correlations fully coincide with 
avalanches of correlated events. This basically means that an 
avalanche occupies a significant fraction of lateral extent of 
the system and cooperative correlated motion over large scales 
does occur. Accordingly, close to the pinning critical point we 
have (£) ~ T 1 / 2 . These scaling relations immediately imply 
(s) ~ T u l z ', so the exponent 7 relating avalanche sizes versus 
durations becomes 



7 = D/z = (d + ai oc )/z, 



(6) 



in the limit v — >• 0. 

This result is to be contrasted with the scaling relation 
obtained in a recent work by Rost et al. [4]. They ana- 
lyzed a regime of relatively high velocities for which the 



3 



length scale £ x is very small as compared with system size 
L. In this regime avalanches are very narrow and one can 
decompose the front motion in independent, spatially local- 
ized, avalanches of forward moves, TJ ava ~ (h/£ iX ) d v. The 
avalanche duration is T ~ w(£)/v(£), where v(£) is the front 
velocity over the region of size I spanned by the avalanche. If 
simultaneous avalanches are narrow and independent events a 
central limit theorem argument gives v(t) ~ g~ d / 2 and this 
leads to (s(T)) - T 1 , with 0] 



7 = ("ioc + d)/(a ioc + d/2), 



(7) 



in d + 1 dimensions. In particular, for d = 1 one has the 
prediction 7 = 4/3 04]]. Interestingly, this argument also 
leads to the scaling relation (£(T)) ~ T s , with an exponent 
S = 1/(q!i oc + d/2) that differs from 1/z, where z = 3 
is the dynamic exponent describing the correlation spread- 
ing of interfacial fluctuations for forced-flow imbibition in 
d = 1 i?, 13 12, 12]. This indicates that the propagation 
of interface correlations is decoupled from the avalanche dy- 
namics. Indeed, as it will be shown below, the scaling the- 
ory leading to Eq. (0 is valid in a velocity regime such that 
the characteristic length scale is negligible as compared with 
the system size £ x <C L. For lower front velocities, when 
£ x becomes comparable with the system size, the exponent 7 
should tend to the value given by Eq. © instead of (|7]). 

Scaling properties in the static limit.- The interface scal- 
ing exponents can be obtained by a scaling theory which is 
expected to be valid in the static (pinned state) limit v — > 0. 
In the pinned state the velocity-dependent term in Eq. ([TJ can- 
cels and the geometric properties of the front can be described 
by the balance between surface tension and capillary disorder, 
which in real space can be written as aV 2 h p {x) + r)(x) ~ 0, 
where h p (x) is the pinned state, and the disorder is delta cor- 
related, (rj(x)ri(x')) — 77Q + A 2 S(x — x'), with a mean value 
(rj(x)) = ijo and variance A 2 — r) 2 ,. Applying a scaling 
transformation, x — > bx and h p —> b a h p , scale-invariance 
holds for a global roughness exponent a — 3/2 for the pinned 
state configuration. Small perturbations of the pinned state Sh 
are assumed to relax towards another of the infinitely many 
pinned configurations according to 



d t {Sh) = aK V 2 {Sh) + K rj(x), 



(8) 



which leads to the exact interface exponents z — 2, a = 3/2, 
and Q!i oc = 1 at the critical point v = for d = 1, These 
exponents can now be replaced in the avalanche scaling re- 
lations (0 and © to obtain tt = t s , n — 2tt — 1, and 
7 = 1, where we have used D = d + a\ oc = 2 in d = 1. In 
the following we compare these scaling results with numeri- 
cal integrations of a phase-field model as one approaches the 
singular point v = 0. 

The phase-field model.- The scaling properties of fluid 
imbibition fronts can be well described by means of a phase- 
field model JUS- A conserved field (f> is used to represent the 
two existing phases, taking the equilibrium values <f> eq = +1 
and <f) eq = —1 in the liquid and air phases, respectively. 
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FIG. 3: (Color online) Average avalanche size and lateral extent ver- 
sus duration in a system of size L — 512 for two velocities in the 
different dynamical regimes discussed in this paper. 



The dynamics of the phase field is controlled by a continuity 
equation based on a time-dependent Ginzburg-Landau model 
with conserved order parameter dt<p = VMV/i where p = 
5 T /5(f) is the chemical potential and the free energy takes the 
form T[4>] = Jdr [(eV</>) 2 /2) - </> 2 /2 + 4 /4 - 7?(r)0] . 
The quenched random field rj(r) > models capillary dis- 
order and favors the liquid (wet) phase, forcing the interface 
to advance at the expense of the air (dry) phase. In our nu- 
merical model we have used a spatially distributed dichotomic 
quenched noise in a two-dimensional system. The locally con- 
served dynamics is described by 



W=VMV[- 



2 V 2 0-7 ? (r)]. 



(9) 



where M is a mobility parameter which we take constant at 
the liquid phase (<j) > 0) and zero at the air phase (<f> < 0), 
and the disorder is Gaussian with a correlator {r]{r)rj(r')) = 
( r l) 2 + (v 2 )3( x ~ x ')- Equation (O is then integrated in the 



"weak" disorder case [11], i.e., when the disorder intensity is 
much smaller than the dimensionless surface tension, in a sys- 
tem size of L = 512 and 25 disorder realizations with e = 1, 
and the forced-flow boundary condition V/i = — vy is im- 
posed at the bottom of the system 01 Oil . All the values for the 
average front velocity have been normalized to a reference 
value vo — 0.002, which corresponds to the highest value 
studied in this paper. 

Numerical results. - Figure|2]shows the avalanche size and 
duration statistics calculated using a threshold c t h = 3. We 
observe that the probability distributions tend to a power-law 
as v is decreased, due to the divergent cutoff [cf. Eq. QJ], 
We estimate the exponents r s ~ 1.54(5) and tt — 1.62(6) 
from the scaling of the data. Despite the smallest velocity we 
were able to reach V = vq/40 = 5 x 10~ 5 , which is still far 
from zero, the scaling region is reasonably good. According 
to our scaling theory both exponents should exactly coincide 
at the critical point, which is consistent with the numerical 
values within the error bars. We also plot a direct estimate 
of the divergent avalanche size cutoff in good agreement with 
Eq. ©. 
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FIG. 4: (Color online) Scaling of the front roughness in the phase- 
field simulations for L — 512 in the low velocity v = wq/40 
limit, (a) Structure factor S(k, t) at different times giving a spec- 
tral roughness exponent of q s = 1.45(6) 11411 . (b) Data collapse of 
the local width according to a super-rough scaling for a — 1.50(3), 
a loc = 1.0, and z = 2.09(6). 



We also find an excellent agreement with our prediction in 
Eq. © for the scaling relation between size and duration of an 
avalanche, (s(T)) ~ T 7 . In Fig.[3]we plot both avalanche size 
and lateral extent vs time for two typical velocities. For the 
lowest velocity we studied v = Urj/40 = 5 X 1CP 5 we estimate 
7 ~ 1.09(7), which is to be compared with 7 = 1 from Eq. (|6} 
in d = 1. This can also be compared with the scaling relations 
between the avalanche exponents r s and tt given by Eq. (0. 
Substituting the numerical values 7 = 1.09 and tt — 1.54 we 
predict r t ~ 1.60 in good agreement with the numerical result 
(cf. Fig. |2). Note that the scaling theory is expected to be exact 
only at the critical point v = 0, which is not actually reached 
with our phase-field model results. However, the singularity 
is strong enough to lead to effective exponents for velocities 
within a critical region v < L~ 2 . 

At variance with Rost et al. Jj], who only monitored 
avalanches in the global velocity time series v(t), here we are 
actually looking at active sites that participate in an avalanche 
and, therefore, we are able to check the validity of the scal- 
ing law {£) ~ T s . The typical lateral extent is predicted to 
scale with avalanche duration with an exponent 8 = l/(a\ oc + 
d/2) = 2/3 in the high velocity regime for d = 1 Jj], in ex- 
cellent agreement with our numerical estimate in Fig. [3] How- 
ever, for low velocities we predict 5 = l/z = 1/2 with the 
dynamic exponent z = 2 in the static limit. A strong proof 
of a distinctive behavior as the front velocity is decreased can 
be readily seen in Fig. [3] Our numerical simulations indicate 
that 7 ->• 1 and (£) ~ T 1 / 2 as v -> 0, pointing out that the 
dynamics is controlled by the static critical point at v = 0. 

On the other hand, we have also estimated the scaling of 
the interfacial fluctuations for states close to the static limit. 
This provides an independent check of the validity of our scal- 
ing relations in Eqs. (01 and (0 connecting avalanche and 
roughness exponents. From the structure factor S(k,t) = 
(\hk(t)\ 2 ) and local width in Fig.|4j we estimate the global, 
local, and spectral roughness exponents 11411 a = 1.50(3), 
aioc = 1, a s = 1-45(6), respectively, and the dynamic ex- 
ponent z = 2.09(6) in excellent agreement with our scaling 
theory for the pinned state. 

Finally, we claim that the existence of a singular behav- 
ior as v — » and the extent of the critical region v < L~ 2 
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TABLE I: Correlated extend and effective scaling exponents from 
phase-field simulations at different velocities in as system of size 
L — 512. The fronts are always super-rough with ot\ oc — 1 fl4Tl . 



explain earlier numerical observations J5I [loll that reported a 
dependence of the critical exponents a(v) and z(v) with the 
velocity in numerical results of forced-flow imbibition in finite 
systems. Table H] summarizes the different interfacial scaling 
exponents we observed for different velocities. We observe 
that as the static limit is approached the dynamic exponent 
z — > 2 and the roughness exponent a — > 3/2, as corresponds 
to the pinned state. 

To conclude, we have studied avalanche dynamics in 
forced-flow imbibition in the pinning limit v — > 0. A scaling 
theory relating the roughness of the front with the avalanche 
dynamics has been developed in excellent agreement with nu- 
merical results. Our scaling analysis is based on the presence 
of long-range correlations due to the divergent characteristic 
length scale at v — > 0. From an experimental point of view, it 
would be of great interest to explore the pinning limit. In the 
experimental setup of a Hele-Shaw cell [6], this limit may be 
achieved by putting the cell at an angle so that gravity plays a 
role. This setup should produce fronts near pinning. Note that 
already in the case of having a correlation length of about a 
30% of the system size, the effect of the critical point should 
show up in a drift of the measured scaling exponents (see 
Table I). Alternatively, fluctuations around the critical point 
v = could also be experimentally tested by setting the cell 
at an angle, so that the front is pinned, and then study how 
the system responds to a small angle variation. In this config- 
uration we expect the front to jump from one pinned state to 
another following a relaxation dynamics described by Eq. ([8j 
driven by avalanches described by Eqs. (O and (|6]). 
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